Phase Field Theory of Heterogeneous Crystal Nucleation 
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A phase field approach is developed to model wetting and heterogeneous crystal nucleation of an 
undercooled pure liquid in contact with a sharp wall. We discuss various choices for the boundary 
condition at the wall and determine the properties of critical nuclei, including their free energy of 
formation and the contact angle as a function of undercooling. We find for particular choices of 
boundary conditions, we may realize either an analog of the classical spherical cap model or decidedly 
non-classical behavior, where the contact angle decreases from its value taken at the melting point 
towards complete wetting at a critical undercooling. 
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Heterogeneous nucleation is not only a phenomenon of 
classic importance in materials science, but also remains 
one of continuously growing interest, due to the emerging 
technological interest in nanopatterning techniques and 
control of related nanoscale processes While solid- 
ification of pure undercooled liquids is initiated by ho- 
mogeneous nucleation (the formation of small crystalline 
fluctuations exceeding a critical size determined by the 
interplay of the driving force of crystallization and the 
interfacial free energy |2(]) the presence of foreign parti- 
cles, container walls, and other heterogeneities typically 
facilitates this process 0|. Despite its vast technological 
importance, heterogeneous nucleation remains poorly un- 
derstood. This deficit stems from the complexity of de- 
scribing the interaction between the foreign matter and 
the solidifying melt. 

Wetting of a foreign wall by fluids/crystals has been 
studied extensively [J] including such phenomena as crit- 
ical wetting and phase transitions at interfaces 0. Vari- 
ous methods have been applied to address these problems 
such as continuum models 0] and atomistic simulations 
0. Despite this inventory, recent studies || address- 
ing heterogeneous crystal nucleation rely almost exclu- 
sively on the classical spherical cap model, which assumes 
mathematically sharp interfaces [9J. Here the wall-liquid 
and wall-solid interactions are characterized by the con- 
tact angle 9 that is determined from the interfacial free 
energies by Young's equation: jwl = Jws + Isl cos(#), 
where subscripts W, S, and L refer to the wall, the solid, 
and the liquid respectively. Such models qualitatively de- 
scribe this system, but lose their applicability Q when 
the size of nuclei is comparable to the interface thick- 
ness (the nanometer range, according to atomistic simu- 
lations @,EJ|). Such nanoscale nuclei are essentially "all 
interface". Recent investigations show that phase 
field theory (PFT.fl^|) can address this issue. Indeed, 
PFT can quantitatively predict the nucleation barrier 
for systems (e.g, hard-sphere, Lennard- Jones, ice-water) 
where the necessary input data are available. We there- 



fore adopt this approach to describe heterogeneous nu- 
cleation. Experimentally, the details of the wall-fluid in- 
teraction are embedded in more directly accessible quan- 
tities, such as the contact angle in equilibrium. It is 
thus desirable to develop a model that describes the wall 
in such phenomenological terms. Along this line, inter- 
action between dendritic growth and wall has recently 
been discussed in while Castro addressed crystal 

nucleation in a specific case of 9 — 90°, obtained by p re- 
scribing "no-flux" boundary condition at the wall [1J|. A 
more general treatment is, however, required. 

In this Letter, we describe how to implement phase 
field methods of heterogeneous nucleation with an arbi- 
trary contact angle. For simplicity, we consider a single 
component system, whose local state is characterized by 
the non-conserved phase field 0(r), defined so that and 
1 correspond to the bulk solid and the liquid phases, re- 
spectively. Following previous work ja, 111- we assume 
that the interaction of the wall with the solidifying sys- 
tem is of sufficiently short range to be characterized by 
a "contact free energy" 7w(</>) that depends only on the 
local state of matter abutting the wall. Accordingly, the 
free energy of the system consists of a surface and a vol- 
umetric contribution 
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Here A is a closed surface bounding volume V of the 
solid-liquid system. At A the system is in contact with 
the wall. The volumetric term is the standard form 
found in thermodynamically consistent formulations of 
PFT 0| . Specifically, the local free energy density, given 
the temperature T, is /(</>) = wTg(^) + [1 -p(<f>)]fs(T) + 
p{.4>) f l(T) , while the "double well" and "interpolation" 
functions have the forms g(<p) — (1/ 4) 2 (1 — (j)) 2 and 
p(4>) — 4> 3 (10 — 150 + 6(f> 2 ). The model parameters can 
be related to the solid-liquid interface free energy, the 
interface thickness S and the melting temperature T m as 
e 2 = 2^ 2 6 1SL S/T m and w = 2 l /H 1S L/{T m 5). 
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The critical fluctuation (nucleus) represents an ex- 
tremum of the free energy. The extremum condition 
SFtot — yields the following equations 



e 2 T(n • V0) =0 on A; 



t 2 T{y 2 (j)) = o 



df 



(2) 
(3) 



Here n is the normal vector pointing away from the wall, 
and Eq. (2) is the boundary condition on A to Eq. 
(3) , which in turn is the differential Euler-Lagrange (EL) 
equation to be satisfied inside the volume V. 

For simplicity, we consider first a semi-infinite system 
in contact with a wall. We label the phase field at the 
wall (f>Q, and in the far field (j>ao. There are three possible 
extrema of the free energy in this case: (i) stable solid in 
contact with the wall (0oo = 0; absolute minimum); (ii) 
metastable liquid in contact with the wall (4>oo = 1 ; local 
minimum); (iii) unstable solid droplet (critical fluctua- 
tion) formed in metastable liquid at the wall (</>oo = 1; 
saddle point). 

To advance further, we must specify the contact free 
energy 7iy(</>), this is, in general, based on the details of 
the wall-fluid interaction. Thus, we now consider two il- 
lustrative choices for our boundary conditions, and relate 
these choices to the equilibrium contact angle. 

Model A: We assume that the wall does not perturb 
the structure of the planar solid-liquid interface. Then 
|V0| can be calculated using the one dimensional version 
of the integral of Eq. (3) (at T = T m ): (e 2 T/2)(V0) 2 = 
f{4>) ~ f{4>oo) — A/((/>), and the normal component of 
the gradient at the wall can be expressed as n • W<j) = 
\V(f>\ cos (8). Combining these expressions we have 

n • = [cos(0)/(2 1 / 2 ( 5)]^ o (l - O ); on A, (4) 

a condition that coincides with for 9 = 90°. The 
respective contact free energy, obtained by integrating 
Eq. (2), reads as yw(4>) -Jwl = ~1sl cos(6>)[2</> 3 -3</> 2 + 
1]. Given the postulated relationship between cos(6>) and 
the interfacial free energies, we find j(<p) — jws — JWL 
at the wall-solid contact and — at the wall-liquid 
contact. We adopt Eq. (4) and the respective jw{4>) m 
the undercooled state. Model A can thus be viewed as a 
phase field implementation of the classical spherical cap 
model similar in spirit to that by Semoroz et al. 

Model B: Alternatively, we may specify (f>o = const.. 
In this case only the EL equation follows from the 
extremum condition. This restriction implies that 
7w(^o) = 7o, independent of whether solid or liquid is in 
contact with the wall. Accordingly, for planar interfaces 
at the melting point, the interfacial free energies can be 
expressed as j W l = h(<f>o,l) + 7o,7ws = h(0, <fo) +70, 



and 7sl = h(0, 1), where 
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(5) 



After some algebra, we find cos(0) = (-/wl— lws)/lLS = 
1 — 6</>Q + 40q. The condition that 4>o — const, sets the 
degree of ordering (or disordering) of the substance abut- 
ting the wall. In other words we see liquid ordering near 
the wall if liquid phase is assumed in the far field, while 
a disordered solid forms at the wall if solid is assumed in 
the far field. Inspection of the integral EL equation indi- 
cates that the metastable liquid solution exists so far as 
A/(0o) > 0, i.e., below a critical undercooling AT C deter- 
mined by the condition Af((f>o) = 0. (As <f>o approaches 
the solid state, AT C converges to 0.) At lower T, there 
is no time-independent solution of this type, instead a 
propagating solidification front emerges that is described 
by the usual equation of motion for the phase field [l5T |. 
This model does not converge to at # = 90°. 

Our two choices of the boundary condition at the wall 
correspond to two distinct physical situations: (a) The 
"classical" case, when liquid ordering is negligible at the 
wall, and (b) a "non-classical" case, where there is an 
imposed order at the wall. This structure is of a spe- 
cific nature, as it corresponds the particular, chosen level 
of ordering as one traverses solid-liquid interface. As 
such, this order is "compatible" with the appearing crys- 
tal structure, and will lower the nucleation barrier to the 
formation of solid. While it is typical for liquids to or- 
der at an abutting wall [7(c), 10(b), 10(d)], such ordering 
may not be compatible with the structure to which the 
liquid crystallizes 0], and a more detailed model would 
be required. Based on these observations, we expect that 
our combined analyses of Models A and B will elucidate 
many of the essential behaviors of physical systems. In 
what follows, we evaluate the properties of heterogeneous 
nuclei in these two limiting cases, and present illustrative 
model simulations for pure Ni hdi. 

The EL equation for the composite system nucleus plus 
undercooled liquid has been solved by the finite element 
method. The initial condition has been created by plac- 
ing the classical sharp interface nucleus into the simu- 
lation window after broadening its interface by a tanh 
function of appropriate width parameter. The simula- 
tion box had the size of 30 nm x 20 nm. The equation of 
motion for dynamic evolution simulations was solved in 
a dimensionless form using the finite difference method 
and parallel computing on a PC cluster of 120 nodes. The 
spatial step was Aa; = 0.2 nm, while noise (as described 
in nas been added to the governing equation. 

Critical fluctuations (nuclei) computed at undercooling 
AT = 35 K as a function of the equilibrium contact angle 
9 are presented in Fig. 1. For comparison, Fig. 2 shows 
the nuclei calculated for a contact angle of 9 = 61.2° 
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FIG. 1: (Color online). Structure of heterogeneous nuclei in 
2D vs. equilibrium contact angle at AT = 35 K in Model A 
(upper row) and B (lower row). Prom left to right 9 = 37.6°, 
72.8°, 107.2°, and 142.3° (c/>o = 0.2, 0.4, 0.6, and 0.8). The 
contour lines stand for (j> = 0.2, 0.4, 0.6, and 0.8. Horizontal 
size is 10 nm. Coloring: red - liquid; yellow - solid. For 
symmetry reasons, only the right half of the nuclei is shown. 
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FIG. 2: (Color online). Structure of nuclei in 2D at three 
undercoolings (AT = 20 K, 40 K, and 90 K) in Models A 
(upper row) and B (central row), at an equilibrium contact 
angle of 9 = 61.2° corresponding to cj>o = 1/3 and AT C = 92.0 
K. The contour lines indicate <f> = 1/6, 2/6, 3/6, 4/6, and 5/6, 
respectively. Horizontal size is 12 nm. Coloring as for Fig. 
1. The free energy of non-classical nuclei (Wa,b) normalized 
by the sharp interface prediction (Wsi) is also shown as a 
function of undercooling (bottom panel): Model A - solid; 
Model B - dashed. 



(^o = 1/3) as a function of AT. While in both models 
the size of heterogeneous nuclei becomes comparable to 
the interface thickness with increasing undercooling, in 
Model B it happens at a far smaller undercooling. It is 
remarkable that while the contact angle is approximately 
constant in Model A, in Model B it varies drastically 



with undercooling in Model B, and tends to 0° complete 
wetting at the critical undercooling. The free energy of 
heterogeneous nuclei in Model A and B are also shown in 
Fig. 2. For Model A, it has been calculated by integrat- 
ing the free energy density difference relative to the bulk 
undercooling liquid and adding the contribution from the 
wall. For Model B, the integrated free energy of the wall- 
liquid system has been subtracted from the free energy 
of the wall-nucleus-liquid system. It is found that the 
free energies of nuclei from Models A and B fall close to 
the values from the sharp interface spherical cap model 
(2D) for small undercoolings where the nuclei are large 
relative to the interface thickness, while lower values are 
obtained at larger undercoolings. In Model B, the nucle- 
ation barrier disappears at AT C . Atomistic simulations 
could test the existence of such a AT C . 

Next, we demonstrate that these models describe wet- 
ting as it is usually understood. For a corner of angle a, 
there exists a critical contact angle 9 C = ir/2 — a/2 below 
which there is no nucleation barrier. For a right angle, or 
square corner, 9 C = 45°. At the melting point we examine 
three cases in a square corner: (a) We select a wetting an- 
gle of 9 = 9 C . An isosceles triangle crystal formed in this 
corner realizes a planar solid-liquid interface. This inter- 
face is stable, since the capillary pressure vanishes for a 
planar interface, (b) For 9 < 9 C , the solid-liquid interface 
is concave, and resulting in a negative capillary pressure 
that increases the melting point (Gibbs-Thomson effect), 
so a crystal triangle posited in the corner must grow, (c) 
For 8 > C , a convex solid- liquid interface develops and a 
crystal triangle placed in the corner melts. Simulations 
performed for rectangular corner display the expected 
behavior for both models (Fig. 3). 

Illustrative simulations demonstrating the power of our 
approach, are presented in Fig. 4. Using Model A, with 




FIG. 3: (Color online). Wetting of rectangular corner by the 
crystalline phase at T = T m in Model A. From up to down 
the contact angle is 40°, 45° (critical), and 50°, respectively. 
Time increases to the right. Snapshots were taken at 5000, 
20,000, 40,000, 60,000 and 80,000 dimensionless time steps. 
(Coloring as for Fig. 1. 100x100 grid.) Similar results were 
obtained for Model B. 
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FIG. 4: (Color online). Phase field simulation of heterogeneous nucleation of Ni crystals on complex surfaces in 3D using Model 
A. From left to right: on stairs, in rectangular grooves, on a checkerboard-modulated surface, and on floating nanoparticles. 
[AT = 300 K. A 400 x 400 x 200 grid corresponding to 80 nm x 80 nm x 40 nm is used. Coloring: solid-liquid interface 
(defined by <j> = 0.5) - red; bulk crystal - black; wall - yellow.] 



9 = 60°, heterogeneous nucleation is modeled on stairs, 
in rectangular grooves, on a checkerboard-modulated sur- 
face, and on floating nanoparticles. Simulations of this 
kind are expected to find applications in nanopatterning 
studies pj. 

We have developed a phase field methodology to de- 
scribe heterogeneous crystal nucleation in undercooled 
liquids at walls characterized by arbitrary contact an- 
gles. Two limiting cases have been addressed: (Model A) 
Nucleation at surfaces where liquid ordering at the wall 
is negligible and (Model B) where the wall-liquid interac- 
tion induces partial crystalline order in the liquid (Model 
B). Using the prescriptions described above, many other 
boundary conditions can be explored. Note that this ap- 
proach can be applied to any system displaying a first 
order phase transition (such as vapor-liquid). Also, this 
approach can be directly extended to existing phase field 
models of alloy and anisotropic polycrystalline systems 
characterized by further fields. The present study thus 
opens up new ways for modeling heterogeneous nucle- 
ation in a broad variety of systems. 
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